#include "CoerciveFieldAlgo.h"

#include "TrajSimul.h"
#include "SwitchingCheck.h"

#include <cmath>
#include <iostream>

CoerciveFieldAlgo::CoerciveFieldAlgo()
{

}

CoerciveFieldAlgo::~CoerciveFieldAlgo()
{

}

double CoerciveFieldAlgo::solve(Clock clock,Output output,Pulse *field,Pulse *current,MultiLayer multiLayer)
{
  Parameters *p=Parameters::getInstance();
  int layerNumber=multiLayer.getSolveLayer();
  int run=p->getInt("nrun");
  if(p->getDouble("temperature")==0)
    run=1;
  double proba=0,nrun=0,nswitch=0;
  double bCritmax=p->getDouble("maxField");
  double bCritmin=p->getDouble("minField");
  double bStep=p->getDouble("fieldStep");
  MultiLayer result;
  bool bmin,bmax;
  field->setLevel(bCritmax);
  MultiLayer tempLayer=TrajSimul().solve(clock,output,field,current,multiLayer);
  multiLayer[layerNumber].setMag(tempLayer[layerNumber].getMag());
  field->setLevel(bCritmin);
  tempLayer=TrajSimul().solve(clock,output,field,current,multiLayer);
  multiLayer[layerNumber].setSwitchedMag(tempLayer[layerNumber].getMag());
  std::cout<<multiLayer[layerNumber].getMag().getZ()<<"=="<<multiLayer[layerNumber].getSwitchedMag().getZ()<<std::endl;
  if(multiLayer[layerNumber].getMag().getZ()==multiLayer[layerNumber].getSwitchedMag().getZ())
    return 0;
  for(nrun=0;nrun<run;nrun++)
    { 
      result=TrajSimul().solve(clock,output,field,current,multiLayer);
      bmin=SwitchingCheck(result[layerNumber]).getCheck();
      if(bmin) nswitch++;
    }
  proba=nswitch/run;
  if(proba>=0.5)
    bmin=true;
  else bmin=false;
      //std::cout<<"switching"<<std::endl;
      //else std::cout<<"not switching"<<std::endl;
  std::cout<<"bmin :"<<proba<<std::endl;
  nswitch=0;
  field->setLevel(bCritmax);
  for(nrun=0;nrun<run;nrun++)
    {
      result=TrajSimul().solve(clock,output,field,current,multiLayer);
      bmax=SwitchingCheck(result[layerNumber]).getCheck();
      if(bmax) nswitch++;
      //std::cout<<"switching"<<std::endl;
      //else std::cout<<"not switching"<<std::endl;
    }
  proba=nswitch/run;
  if(proba>=0.5)
    bmax=true;
  else bmax=false;
 std::cout<<"imin :"<<proba<<std::endl;
 if(!bmin||bmax)
    {
      std::cout<<"bmax ou !bmin"<<std::endl;
      return 0;
    }
 while(fabs(bCritmax-bCritmin)>bStep)
    {
      bool btemp;
      double bTemp=(bCritmax+bCritmin)/2;
      field->setLevel(bTemp);
      nswitch=0;
      for(nrun=0;nrun<run;nrun++)
	{
	  result=TrajSimul().solve(clock,output,field,current,multiLayer);
	  btemp=SwitchingCheck(result[layerNumber]).getCheck();
	  if(btemp) nswitch++;
	}
      proba=nswitch/run;
      if(proba>=0.5)
	btemp=true;
      else btemp=false;
      if(btemp&&bmin)
	{	
	  bCritmin=bTemp;
	  bmin=true;
	}
      if(!btemp&&!bmax)
	{	
	  bCritmax=bTemp;
	  bmax=false;
	}
    }
  return bCritmax;
}
